##########################################################################################

library(Seurat)
library(data.table)
library(optparse)
library(ggplot2)
library(dplyr)
library(RColorBrewer)

##########################################################################################

option_list <- list(
	make_option(c("--sample_list_file"), type = "character"),
    make_option(c("--singleCell_sample_file"), type = "character"),
    make_option(c("--single_cell_scissor_file"), type = "character"),
    make_option(c("--single_cell_file"), type = "character"),
    make_option(c("--gene"), type = "character"),
    make_option(c("--out_path"), type = "character")
)

if(1!=1){

    gene <- "GKN1_GKN2"
    work_dir <- "~/20220915_gastric_multiple/dna_combinePublic/"
    singleCell_sample_file <- "~/20220915_gastric_multiple/dna_combinePublic/config/singleCell_Sample.useThree.list"
    single_cell_file <- paste0(work_dir,"/public_ref/singleCell/njmu/epiall_nor_PCA_50_RE0.5.Rdata")
    single_cell_scissor_file <- paste0(work_dir,"/images/singleCell_MUC6/Scissor_STAD_MUC6_mutation.IM.CellRate.all.RData")
    out_path <- paste0("~/20220915_gastric_multiple/dna_combinePublic/images/singleCellRatio/")
    sample_list_file <- "~/20220915_gastric_multiple/dna_combinePublic/config/tumor_normal.class.MSS_MSI.list"

}

##########################################################################################

parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

sample_list_file <- opt$sample_list_file
out_path <- opt$out_path
singleCell_sample_file <- opt$singleCell_sample_file
gene <- opt$gene
single_cell_file <- opt$single_cell_file
single_cell_scissor_file <- opt$single_cell_scissor_file

dir.create(out_path , recursive = T)

##########################################################################################

info_singlecell <- data.frame(fread(singleCell_sample_file))
info <- data.frame(fread(sample_list_file))

## 总的单细胞
sc_dataset <- load(single_cell_file, verbose = F)
sc_dataset_all <- epiall_nor_PCA_50_RE0.5
##Idents函数定义要取的类别
Idents(sc_dataset_all) <- "sample"   

## scissor分类过
sc_dataset_scissor <- load(single_cell_scissor_file, verbose = F)
sc_dataset_scissor <- sc_dataset

##########################################################################################
## scissor阳性细胞
muc6_pit <- names(which(sc_dataset_scissor$scissor=="1" & sc_dataset_scissor$celltype=="Pit"))
other_pit <- names(which(sc_dataset_scissor$scissor!="1" & sc_dataset_scissor$celltype=="Pit"))

##########################################################################################
## 提取最后纳入分析的样本
## 9个

info_singlecell <- subset( info_singlecell , singlecell_ID != "" )
info_singlecell <- info_singlecell[info_singlecell$ID %in% info$ID ,]

##########################################################################################
## 计算细胞比例
## 合并在一起
result_tmp <- c()

for( class in c("IM") ){

	print(class)

	## 提取特定病理类型的样本
	sc_dataset <- subset(sc_dataset_all , idents=c(class))
	Idents(sc_dataset) <- "celltype"   
	## 只看Pit细胞
	sc_dataset <- subset(sc_dataset , idents=c("Pit"))

	## 去除MUC6突变样本
	info_singlecell_use <- subset( info_singlecell , Class == class & ID != "JZGC00732" )

	# 计算每个细胞中基因的表达量
	expression_matrix <- as.matrix(GetAssayData(sc_dataset, assay = "RNA"))
	## 中位表达
	he_t1 <- as.numeric(quantile(expression_matrix["GKN1",])["50%"])
	he_t2 <- as.numeric(quantile(expression_matrix["GKN2",])["50%"])

	## 计算每个人中高表达基因的比例

	for( sample in info_singlecell_use$singlecell_ID ){
		index <- grep( sample , colnames(expression_matrix) )

		tmp_expression <- expression_matrix[c("GKN1" , "GKN2"),index]
		high_ratio <- length(which(tmp_expression["GKN1",] > he_t1 & tmp_expression["GKN2",] > he_t2))/ncol(tmp_expression)

		tmp <- data.frame( Sample = sample , High_Ratio = high_ratio )
		result_tmp <- rbind(result_tmp , tmp)
	}

	## scissor鉴定的突变pit
	tmp_expression <- expression_matrix[c("GKN1" , "GKN2"),muc6_pit]
	high_ratio <- length(which(tmp_expression["GKN1",] > he_t1 & tmp_expression["GKN2",] > he_t2))/ncol(tmp_expression)
	tmp <- data.frame( Sample = "JZ732P\n(MUC6 Mut)" , High_Ratio = high_ratio )
	result_tmp <- rbind(result_tmp , tmp)

	## scissor鉴定的其它pit
	tmp_expression <- expression_matrix[c("GKN1" , "GKN2"),other_pit]
	high_ratio <- length(which(tmp_expression["GKN1",] > he_t1 & tmp_expression["GKN2",] > he_t2))/ncol(tmp_expression)
	tmp <- data.frame( Sample = "JZ732P\n(Other)" , High_Ratio = high_ratio )
	result_tmp <- rbind(result_tmp , tmp)

}


result_tmp_lowratio <- data.frame( Sample = result_tmp$Sample  , High_Ratio = 1 - result_tmp$High_Ratio)
result_tmp_lowratio$Type <- "Low"
result_tmp$Type <- "High"

result_final <- rbind( result_tmp , result_tmp_lowratio )
colnames(result_final)[2] <- "Ratio"

## 样本排序
sample_order <- result_tmp$Sample[order(result_tmp$High_Ratio , decreasing=T)]
result_final$Sample <- factor( result_final$Sample , levels = sample_order , order = T )
result_final$value_text <- paste0( round(result_final$Ratio , 2) * 100 , "%") 

col <- c(
  brewer.pal(9,"YlGnBu")[6],
  rgb(255,0,0,alpha=255,maxColorValue=255)
  )
names(col) <- c("Low" , "High")

p1 <- ggplot(result_final,aes(x=Sample,y=Ratio,fill=factor(Type))) +
	geom_bar(stat="identity") +
	ylab(paste0( "High expression cells proportion of\n" , gene)) +
	geom_text(aes(label=value_text) , position=position_stack(vjust = 0.5) , size=3.5 , color="black")+
	xlab(NULL) +
	theme_bw() +
  	theme(panel.background = element_blank(),#设置背影为白色#清除网格线
        legend.position ='top',
        legend.title = element_blank() ,
        panel.grid.major=element_line(colour=NA),
        legend.text = element_text(size = 8,color="black",face='bold'),
        axis.text.x = element_text(size = 10,color="black",face='bold'),
        axis.text.y = element_text(size = 8,color="black",face='bold'),
        axis.title.x = element_text(size = 10,color="black",face='bold'),
        axis.title.y = element_text(size = 10,color="black",face='bold'),
        strip.text.x = element_text(size = 12,color="black",face='bold',angle=45,vjust=0.5,hjust=0.5),
        axis.line = element_line(size = 0.5))  +
	scale_fill_manual(values=c(col))


out_name <- paste0( out_path , "/IM_Pit_HighExpressionCellRatio." , gene , ".pdf"  )
ggsave(file=out_name,plot=p1,width=9,height=4)


##########################################################################################
## 其它样本的合起来
result_tmp <- c()

for( class in c("IM") ){

	print(class)

	## 提取特定病理类型的样本
	sc_dataset <- subset(sc_dataset_all , idents=c(class))
	Idents(sc_dataset) <- "celltype"   
	## 只看Pit细胞
	sc_dataset <- subset(sc_dataset , idents=c("Pit"))

	info_singlecell_use <- subset( info_singlecell , Class == class )

	# 计算每个细胞中基因的表达量
	expression_matrix <- as.matrix(GetAssayData(sc_dataset, assay = "RNA"))

	index1 <- grep( info_singlecell_use$singlecell_ID[1] , colnames(expression_matrix) )
	index2 <- grep( info_singlecell_use$singlecell_ID[2] , colnames(expression_matrix) )
	index3 <- grep( info_singlecell_use$singlecell_ID[3] , colnames(expression_matrix) )
	index4 <- grep( info_singlecell_use$singlecell_ID[4] , colnames(expression_matrix) )

	expression_matrix <- expression_matrix[,c(index1,index2,index3,index4)]

	## 中位表达
	he_t1 <- as.numeric(quantile(expression_matrix["GKN1",])["50%"])
	he_t2 <- as.numeric(quantile(expression_matrix["GKN2",])["50%"])

	## scissor鉴定的突变pit
	tmp_expression <- expression_matrix[c("GKN1" , "GKN2"),muc6_pit]
	high_ratio <- length(which(tmp_expression["GKN1",] > he_t1 & tmp_expression["GKN2",] > he_t2))/ncol(tmp_expression)
	tmp <- data.frame( Sample = "MUC6\nMut" , High_Ratio = high_ratio )
	result_tmp <- rbind(result_tmp , tmp)

	## scissor鉴定的其它pit
	tmp_expression <- expression_matrix[c("GKN1" , "GKN2"),colnames(expression_matrix)[!(colnames(expression_matrix) %in% muc6_pit)]]
	high_ratio <- length(which(tmp_expression["GKN1",] > he_t1 & tmp_expression["GKN2",] > he_t2))/ncol(tmp_expression)
	tmp <- data.frame( Sample = "Other" , High_Ratio = high_ratio )
	result_tmp <- rbind(result_tmp , tmp)
	
}


result_tmp_lowratio <- data.frame( Sample = result_tmp$Sample  , High_Ratio = 1 - result_tmp$High_Ratio)
result_tmp_lowratio$Type <- "Low"
result_tmp$Type <- "High"

result_final <- rbind( result_tmp , result_tmp_lowratio )
colnames(result_final)[2] <- "Ratio"

p <- fisher.test(round(matrix(c(result_final$Ratio),ncol=2) * 100))$p.value
trans <- function(num){
    up <- floor(log10(num))
    down <- round(num*10^(-up),2)
    text <- paste("p == ",down," %*% 10","^",up)
    return(text)
}
if( p < 0.01 ){
    p_text <- trans(p)
}else{
    p_text <- paste0( "p == " , round(as.numeric(p) , 3) ) 
}

## 样本排序
result_final$value_text <- paste0( round(result_final$Ratio , 2) * 100 , "%") 
result_final$value_text <- paste0( round(result_final$Ratio , 2) * 100 , "%") 

col <- c(
    rgb(red=179,green=34,blue=35,alpha=255,max=255), 
    rgb(red=2,green=100,blue=190,alpha=255,max=255) 
    )

names(col) <- c("High" , "Low" )

p1 <- ggplot(result_final,aes(x=Sample,y=Ratio,fill=factor(Type))) +
	geom_bar(stat="identity") +
	ylab(paste0( "Proportion of cells with high expression \nGKN1 and GKN2 in pit cells")) +
	geom_text(aes(label=value_text) , position=position_stack(vjust = 0.5) , size=4 , color="white")+
	geom_text(aes(label=p_text , y = 1.05 , x = 1.5),parse = TRUE,size=5 , color = "black") +
	xlab(NULL) +
	theme_bw() +
  	theme(panel.background = element_blank(),#设置背影为白色#清除网格线
        legend.position ='right',
        legend.title = element_blank() ,
        panel.grid.major=element_line(colour=NA),
        legend.text = element_text(size = 12,color="black",face='bold'),
        axis.text.x = element_text(size = 12,color="black",face='bold'),
        axis.text.y = element_text(size = 10,color="black",face='bold'),
        axis.title.x = element_text(size = 14,color="black",face='bold'),
        axis.title.y = element_text(size = 14,color="black",face='bold'),
        strip.text.x = element_text(size = 12,color="black",face='bold'),
        axis.ticks.length = unit(0.2, "cm") ,
        axis.line = element_line(size = 0.5)
        )  +
	scale_fill_manual(values=c(col))

print(result_final)
out_name <- paste0( out_path , "/IM_Pit_HighExpressionCellRatio." , gene , ".combine.pdf"  )
ggsave(file=out_name,plot=p1,width=3.5,height=4.5)
